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The stability of stratified, rotating systems 
and the generation of vorticity in the Sun 

Steven A. Balbus 1,2 ' 3 , Emmanuel Schaan 1,3 
ABSTRACT 

We examine the linear behavior of three-dimensional Lagrangian displace- 
ments in a stratified, shearing background. The isentropic and iso-rotation sur- 
faces of the equilibrium flow are assumed to be axisymmetric, but otherwise fully 
two-dimensional. Three-dimensional magnetic fields are included in the perturba- 
tion equations; however the equilibrium is assumed to be well-described by purely 
hydrodynamic forces. The model, in principle very general, is used to study the 
behavior of fluid displacements in an environment resembling the solar convection 
zone. Some very suggestive results emerge. All but high-latitude displacements 
align themselves with the observed surfaces of constant angular velocity. The 
tendency for the angular velocity to remain constant with depth in the bulk of 
the convective zone, together with other critical features of the rotation profile, 
emerge from little more than a visual inspection of the governing equation. In the 
absence of a background axial angular velocity gradient, displacements exhibit 
no poleward bias, suggesting that solar convection "plays-off " of prexisting shear 
rather than creates it. We argue that baroclinic vorticity of precisely the right 
order is generated at the radiative/convective zone boundary due to centrifugal 
distortion of equipotential surfaces that is not precisely followed by isothermal 
surfaces. If so, many features of the Sun's internal rotation become more clear, 
including: i) the general appearance of the tachocline; ii) the extension of differ- 
ential rotation well into the radiative zone; iii) the abrupt change of morphology 
of convective zone isorotation surfaces; and iv) the inability of current numeri- 
cal simulations to reproduce the solar rotation profile without imposed entropy 
boundary conditions. 
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1. Introduction 

In its most general form, the dynamical state of the interior of a star is one of differential 
rotation and entropy stratification. If isobaric and isochoric surfaces do not coincide, the 
angular velocity need not be constant on cylinders. A noteable example is the Sun, for which 
helioseismology studies have fashioned a remarkably detailed and rich portrait. Where the 
Sun is stably stratified in entropy, in the bulk of the radiative zone, it tends not to be 
differentially rotating. However, surrounding the radius of vanishing entropy gradient, in 
both the convective and radiative layers, there is significant differential rotation generally 
dominated by the radial component of the angular velocity gradient. Higher in the convection 
zone, the rotation contours show an abrupt change in morphology, with the sudden emergence 
of a distinctly conical pattern (e.g. Miesch & Toomre 2009). Finally, approaching the Sun's 
surface, there is once again an abrupt shift in contour morphology, apparently associated 
with the onset of high-velocity convection. 

In previous work (Balbus, Latter, & Weiss 2012 [BLW] and references therein), it has 
been shown that the pattern of coaxial cones in the bulk of the solar convective zone (SCZ) 
can be understood as an elementary solution of the vorticity equation (in the limit of thermal 
wind balance) under certain well-posed assumptions. One of these assumptions is that a small 
angular entropy gradient is present, for without this there can be no axial component of the 
angular velocity gradient. Where does this all-important entropy gradient come from? Is 
it, as is often argued, an ineluctable consequence of convection and the Coriolis force, or 
is something more — or something else — involved? For that matter, is the entropy gradient 
more or less fundamental than the concomitant angular velocity gradient? 

To address these questions, we begin with a very general study of the linear behavior 
of three-dimensional fluid displacements in a shearing and stratified background medium. 
The background angular velocity and entropy profiles may depend upon both poloidal co- 
ordinates. It is demonstrated that for a medium in uniform rotation, the most unstable 
displacements do not deflect from spherically radial paths, despite the presence of Coriolis 
forces. When a axial component of the angular velocity gradient is already present however, 
it is shown that there is a significant polar deflection of higher entropy fluid elements. More 
precisely, the sense of this deflection is poleward for outward-moving displacements if the 
axial angular velocity gradient is negative (as in the Sun), and equatorial if this gradient is 
positive. Thus, a axial angular angular velocity gradient is self-reinforcing, and may thus be 
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reshaped, in a convective fluid. Angular velocity gradients in cylindrical radius are much less 
effective in this regard: they have no first order effect on convective displacements. Because 
a background axial angular velocity gradient requires a vorticity source, this has far reaching 
consequences, and we use this finding as a thin edge of wedge to pry further into the origins 
of the Sun's baroclinic differential rotation. We argue in particular that it is likely that the 
rotation pattern of the SCZ has emerged by responding to a preexisting angular entropy 
gradient, rather than generating such a gradient internally. This is entirely consistent with 
the experience of numerical simulations, in which rotation on cylinders stubbornly persists 
unless latitudinal entropy boundary conditions are present, in which case solar-like profiles 
emerge relatively easily (e.g. Miesch, Brun, Toomre 2006). Indeed, if the conclusions of this 
paper are well-founded, the direct imposition of such boundary conditions is the "correct" 
procedure! 

In the second part of this paper, we put forth the case that vorticity generation is 
all but inevitable near the outer edge of the radiative zone where the entropy gradient 
vanishes. The combination of diffusive heating and centrifugal distortion of equipotential 
surfaces is incompatible with radiative and dynamical equilibrium in a uniformly rotating 
medium. The classic remedy of introducing a tiny amount of meridional circulation (e.g., 
Schwarzschild 1958) breaks down at a surface of zero entropy gradient. Instead, radiative 
equilibrium is re-established with very slightly different isothermal and isochoric surfaces. If, 
as one would expect from radiative considerations, the isotherms are more spherical than the 
isochoric surfaces, an incipient tachocline is generated, bearing many of the features observed 
of the true solar tachocline: a negative axial gradient of the angular velocity everywhere, 
a dominant (spherical) radial component of this gradient, and an increasingly dominant 
cylindrical disposition of the isorotation contours toward the equator. 

The generation of vorticity at a level stemming from the centrifugal distortion of the 
equipotential surfaces has not been hitherto viewed as an important component of the solar 
differential rotation profile. However, not only is the centrifugal distortion of precisely the 
correct order-of-magnitude for this problem, we are argue that it is the principal causal 
agent. If this is correct, numerical simulations whose goal is to reproduce the Sun's internal 
rotation accurately from first principles will ultimately have to accomodate this 1 : 10 5 effect. 

The organization of the paper is as follows. In §2, we present the governing equations 
for three-dimensional fluid displacements in an axisymmetric but otherwise fully general 
entropy-stratified, shearing background. This is an interesting gasdynamical problem in its 
own right, and particularly relevant for the sun. General solutions are presented in §3, but 
we focus on the most rapidly growing modes, for which a simple analysis is possible. The 
solution shows explicitly the relationship between shear, Coriolis forces, and the deflection of 
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convective trajectories. In §4, we integrate our findings with standard solar models, arguing 
that the seed angular entropy gradient is a result of centrifugal distortion of equipotential 
surfaces in the radiative zone together with the disappearance of the entropy gradient at the 
SCZ boundary Finally, §5 summarizes our results. 



2. Linear convection theory: fundamental equations 

2.1. Equilibrium state 

Throughout this paper, we use standard cylindrical coordinates (with R the radial 
distance from the rotation axis, the azimuthal angle, and z the distance along the rotation 
axis) and standard spherical coordinates (with r the radial distance from the origin, 9 the 
colatitude angle from the z axis, and <fi the azimuthal angle). Unit vectors will be denoted 
by an appropriately subscripted e. 

The unperturbed background state is one of time-steady hydrostatic equilibrium, 

Ml 2 e H =-VP + V$ (1) 
P 

Here, Q is the angular velocity, p the mass density, P the gas pressure, and $ the gravitational 
potential. The magnetic field B is assumed to be weak and may be ignored in the equilibrium 
state, but large wavenumber perturbations could in principle be significantly influenced by 
the field. We therefore will retain the magnetic field when analyzing small disturbances. 
Equation ([TJ) is quite general for hydrostatic stars or disks, but for SCZ applications it is an 
excellent approximation to take $ = —GM Q /r (G is the Newtonian constant and M & is one 
solar mass), and to ignore the centrifugal force in the equilibrium state. 

The centrifugal term can of course never be ignored in the vorticity equation (the 
component of the curl of equation [I]), which provides a measure of the departure of isobaric 
and isochoric surfaces: 

dQ 2 1 fdpdP dpdP s 

which may also be written 

dz ^p \ dz dR dR dz 

where 7 is the adiabatic index, and a = In Pp -7 is proportional to the specific entropy. The 
right side of (J3J) may be written in spherical coordinates as 

R dQ 2 1 (dadP dadP\ _ g_d_o_ approx i mat i on ) (4) 

dz r^p \ dr d9 d6 dr J d9 



R dz p 2 \dRdz dzdR/' U) 



R dVl 2 _ 1 fdadP da dP^, (0) 
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where g = —(l/p)(dP/dr) is an excellent approximation to the gravitational acceleration, 
and the final approximate equality assumes that the radial component of the entropy gradient 
does not exceed the latitudinal component by many orders of magnitude. Equation (BJ, often 
referred to as the thermal wind equation (e.g. Pedlosky 1987), appears to be well satisfied 
throughout much of the SCZ. 



To understand more fully the complex dynamics of the SCZ, we analyze here a much 
simpler proxy system: the local linear behavior of three-dimensional linear disturbances 
in weakly magnetized, stratified, differentially rotating, two-dimensional background flows. 
Note that although the equilibrium is axisymmetric, the disturbances are fully three-dimensional. 
Our focus is the temporal behavior of the fluid displacements embedded in such a medium. 

Nonlinear convection generally involves coherent, extended structures. A local WKB 
treatment cannot hope to capture fully this element of the problem. Instead, by affording 
some insight as to how uniformly rotating surfaces arise and host convective displacements, 
the local theory can suggest the origin of structure on larger scales. Such structure is able 
to survive despite the presence of shear. 

The fundamental dynamical equation of motion for the linear perturbations is 



As noted, while the magnetic field B is assumed to play no role in the equilibrium state, 
it can still be important for the evolution of large wavenumber perturbations, and will be 
retained. We wish to explore the nonaxisymmetric behavior of local disturbances. Because 
the background flow is in a state of differential rotation, we work in a locally shearing 
Lagrangian coordinate system, and the linear perturbations are ultimately to be expressed 
in terms of the Lagrangian fluid element displacement £(R, <f), z, t). 

Begin by taking standard Eulerian perturbations (5v, 5P, etc.) of the usual fluid equa- 
tions. The linearized equation of motion is 



2.2. 



Inertial terms 





D 



Dt 



5v + (5v-'V)v 




B5B 

Air 



(B-V)5B 

4np 



(6) 



where 




Dt dt d(p 



(7) 
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is the Lagrangian time derivative associated with the unperturbed flow. The magnetic pres- 
sure buoyancy term and the magnetic tension term involving the gradient of the background 
magnetic field have been dropped under the assumption that the field is weak. We make 
the standard WKB assumption that the product of the perturbation wavenumber with any 
background scale height is large. 

Consider the first two terms of equation (JSJ), 

D dSv dSv 

—Sv + (SvV)v = ^- + SI— + (5vX7)(me^ (8) 

Dt at o<p 

The R and <fi components of these terms are, respectively, 



D 
Dt 



Sv + (Sv-'V)v 



-j=^$ v + (Sv-^V)v 



£6vr 
Dt 

.2 



2Q Sv, 



where 



2Q 



n 



d(RQ) 
OR 



The z component is simply 



—Sv + (Sv-'V)v 

Dt K ' 



D5v z 
Dt 



(9) 
(10) 

(11) 
(12) 



The relationship between £ and Sv is given by 
For j equal to R or z, we find Svj = D^j/Dt, while the component of (fT3"j) gives 



Dt 



Sv^ + R{£-V)Q. 



(13) 



(14) 



Equations (jHJ) and ( fTUI) may now be expressed in terms of £r and becoming respectively 

D 



Dt 



Sv + (Sv-'V)v 



e<t>- 



D 
Dt 



Sv + (Sv-'V)v 



U + 25^, 



(15) 
(16) 
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where the dot notation indicates the Lagrangian derivative D/Dt. Putting the last two 
equations together with the z equation of motion leads to the vector equation 

^-Sv + (5vV)v = £ + 2Q x £ + e R R(£-V)n 2 (17) 

The second term on the right is obviously the Coriolis term, and the final term is the "residual 
centrifugal force:" the difference between the centrifugal force in the rotating frame and the 
forces maintaing the differential rotation. Notice the appearance of Vfi 2 , as opposed to an 
angular momentum gradient, as part of the inertial forces. 



2.3. Linear perturbations: magnetic induction and entropy constraints 

Next, recall the relationship between SB and the displacement which follows from 
the integrated form of the induction equation, 

SB = Vx(£xB) (18) 

We work in the Boussinesq limit, 

V-£ = 0, (19) 

so that equation ffl8|) becomes 

SB = (B-V)t (20) 

We have dropped the term (£-V)£? under the WKB assumption that the displacements are 
rapidly varying in space. 

Finally, for adiabatic perturbations, 

7- = £■ Vor, (21) 

p 

since in the Boussinesq limit, the relative pressure perturbation SP/P is small compared 
with the relative density perturbation Sp/p. 



2.4. Comoving coordinates and wavenumbers 

2.4-1- Time dependence of Eulerian wavenumbers 

The final step is to transform to spatial Lagrangian coordinates comoving with the 
unperturbed flow, the "primed" coordinate system. This is accomplished by the transfor- 
mation 

B! = R, 4>' = 4>- Qt, z = z, t' = t, (22) 
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where of course fHs a function of R and z. Use of the Lagrangian derivative D/Dt (= d/dt') 
has already effected the transformation of the partial time derivative, and the two altered 
poloidal spatial derivatives are 

°- = -?--t™°- (23) 
OR OR' dRd^ 1 ; 

d d _ dn d 

~d~z ~ ~dz~> _t 9J90 7 ' ^ ' 

More compactly, 

V = V'-(tVfi)^, (25) 

a relation that holds for all three components of the gradient. In the WKB limit, all distur- 
bances in Lagrangian comoving coordinates have the spatial dependence 

exp [i (k' R R' + m<f>' + k' z z')\ (26) 

in which all components of the k' wave vector are constants. This means that the Eulerian 
spatial derivatives of R and z are replaced locally (and respectively) by ik R (t) and ik z (t), 
where 

k R (t) = k' R -mt—, (27) 
dQ 

k z (t) = k' z - mt— (28) 
Henceforth, the time dependence of kn and k z will be understood with 

■ dVL 

k R = -m— , k z = -m—. (29) 

Note that for initially purely azimuthal e tm ^' disturbances, the poloidal wavenumber com- 
ponents comprise a two-dimensional vector parallel to — Vf2. 



2.4-2. Lagrangian displacements and isorotation surfaces 

At sufficiently large times t (or all times if k' R — k' z — 0), 

k R dn/dR 

T z ^ ~dn/dz~ m 

and the wave vector becomes increasingly axisymmetric as the poloidal components grow. If 
the three components of the displacement £ are of comparable magnitude, then the condition 
k • £ = becomes, at sufficiently large t, 

£-Vfi = 0, (31) 
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whence 

i-Vtt = 0. (32) 

In other words, independently of the details of the dynamics, the velocity vector of a dis- 
turbed fluid element must eventually lie in a surface of constant f2. There is nothing "solar" 
about this argument, it relies entirely on the kinematics of differential rotation and mass 
conservation in an incompressible fluid. But it is tempting to apply this to the Sun, since 
the fluid elements in question would then be convecting heat and eliminating excess entropy 
gradients within surfaces of constant Q, and the confluence of these surfaces with constant 
residual entropy surfaces becomes less mysterious. There is a further benefit to such an ap- 
proach: the preponderance of the constant Q surfaces are observed to be quasi-radial spokes, 
as seen in merdional cross section. They are customarily if crudely described as cones of 
constant 9. If the radially convecting elements are compelled to follow paths of constant Vt, 
then VL ~ Vt(9) is hardly mysterious. (Precisely the same vanishing-divergence reasoning also 
leads to the conclusion that the poloidal components of the magnetic field vector should lie 
in constant Q surfaces.) 

Is this kinematical argument for the alignment of constant residual entropy and angular 
velocity surfaces correct? This depends upon whether there is sufficient time for the shear 
to shape the wavelet before coherence is lost. A potential difficulty is that, as the differential 
rotation is not large, this may be a rather long time interval depending upon the initial 
poloidal wavenumber components. Long-lived, coherent, nonlinear and nonlocal structures 
are seen in fully-developed convection, and these will in fact tend to lie in constant Q surfaces. 
In some sense this backs our approach, but at the same time it puts the cart before the horse: 
it is just this rotation-entropy link we would like to understand. We suggest that a more 
rapid linear dynamical explanation is also available, whereby modes with very small initial 
poloidal wavenumbers and near radial displacements are preferred for rapid development. 
This is discussed in more detail in §2.6 and §3. 



2.4-3. The constancy of k • B 

The component of the equilibrium magnetic field is not independent of time, but 
satisfies the induction equation 

^ = R(B.Vn (33) 

or 

B (j> = B'+tR(B-V)a (34) 
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Note, however, that the magnetic tension k • B — k' • B f , where B' is the magnetic field at 
t — 0, is independent of time. Introducing the Alfven velocity 

VA s vfe' (35) 

the quantity (k • va) 2 , the local magnetic tension force per unit mass, may be regarded as 
a locally constant parameter in the equations of motion. 



2.5. Final Dynamical Equations 

We may now assemble the three fundamental dynamical equations of motion: 

in - 2^ + R(t-V)tf - + ik R (— + *^**) + (k • v A fi R = 0, (36) 

oR 7p \ p Amp J 

& + 2% + ^(— + + (fc • v A )% = 0, (37) 



R \ p Airp 

_ d P(^ + ^ UP + B_6B^ + . = Q 

oz ^p \ p Airp J 

The three dynamical equation may be combined into a single vectorial equation 

£ + 2ftx£ + e^(£-V)ft 2 - — ^Vff + (fc • VA ) 2 £ + — ( 5P + ^t^-) = (39) 

IP P \ 4vr 7 

with the understanding that fc in the final term is time dependent, 

k(t) = fc(0) - mtVQ, (40) 

and 

k(t) • £ = 0. (41) 



2.6. Self-consistent radial convection 



It is a curious and significant fact that in the bulk of the convective zone the Sun tends to 
eliminate every extraneous nonconvective term in the linear equation (1391) . either by the term 
vanishing identically or by its cancellation with another nonconvective term. The dominant 
dynamics is due almost entirely to the underlying radial forcing from the unstable entropy 
gradient, even in the presence of rotation. 
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Recall the concept of residual entropy introduced by Balbus et al. (2009): the entropy 
<r(r, 9) is written as the sum of a function depending only upon spherical radius, a r (r), and a 
residual term, a'(r,6). Physically, a r represents the underlying convection-driving unstable 
radial entropy profile, and a' is the ^-dependent modification that results as a consequence 
of rotation plus convection. In numerical simulations, this breakdown has an operational 
significance: a r is externally imposed, and a' is, in essence, the computed reponse (Miesch et 
al. 2006)0. Only a' is relevant to the thermal wind equation (j2J), since the entropy appears 
exclusively in the form of da/ 89. 

We write a = a r + a' and let us assume that P is a function of r only. Then with 
ge r = — (l/p)VP, the usual squared Brunt- Vaisala frequency N 2 is 

„ r9 q da q d(a r + a') - q da' . . 

N = 'IT = " = N r + - o-- 42 

7 ar 7 or 7 or 

In the SCZ, N 2 < 0. Dropping the magnetic terms, as they appear to be genuinely tiny, 
equation (1591 may then be written 

VP ikS P 

£ + N%e r + 2fix£ + e R R(£-V)n 2 £-Va' + = (43) 

IP P 

Recall that the bulk of the SCZ is characterized by an angular velocity that is insensitive 
to depth; roughly speaking, Q ~ Q(9). Moreover, surfaces of constant fl and a' coincide well 
(Miesch et al. 2006; Balbus et al. 2009). Under these circumstances, for the dominant radial 
e r convective displacements, all terms in equation ( fj^ ) vanish, cancel hydro statically, or are 
otherwise negligible, save the first and second. In particular, the Coriolis deflection (third 
term in from the left) is balanced by the azimuthal pressure gradient (final term on the left 
side), and the two £ ■ V terms are are intrinsically small. The radially moving disturbances 
disturbances are characterized by poloidal wavenumber components very small compared 
with m/R. In other words, within the context of simple linear theory in a uniformly rotating 
sphere, we have a plausible beginning for understanding why the Sun's gross pattern of 
differential rotation is the way it is: the most efficient way to convect heat outwards is 
by maintaining radial convection, which is however permitted only to the extent that the 
dynamical forces of differential rotation allow it. With Q = fl(9), and fluid motions embedded 
within coinciding Q and a' surfaces, the dynamics of radial covection in a shearing system is 
self-consistent. In fact, the data show that in the bulk of the SCZ, constant Q surfaces are 
slightly more axial than constant 9 surfaces. In the next section, we will see that poleward 



1 In practice, the computed a' may acquire a purely spherical contribution as well, but this is easily 
removed by subtracting off the mean. 
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trajectory deviations emerge from the solutions of equation (j43p when dVt/dz < 0. We shall 
argue, moreover, that these axial departures from radial trajectories furnish an important 
clue to the origin of the Sun's vorticity. 



2.7. Reduction to Two Coupled Equations 



The equation of mass conservation V-£ = may be written 

R 



U = ( k R ^ R + k ziz) ■ 

m 



From this it follows 
R 



m 



k R U + Ki z + v)fi, 



R 



rn 



(44) 



k R U + k z U+2R(i.V)n. (45) 



Next, from equation (1381) . 

fSP B • 8B 

\ p Airp 



1 



2t dP(^V)a 



(46) 



dz 7p 

Substituting equations ( 1441) - ( 1461) into ( 1361) and ( 1371) and simplifying produces the equations 

k R f- z „ \ 2QR f. . . : \ DP 



^ + (fc • ^fe - (6 + (fe • ^a) 2 6 



k R ^ R + k z U + (^VV = 0, (47) 

m v / p7 



m 



lfl + (*:-«A) 2 e 



+ 



m Rk, 



& + (fe • »a) 2 & 



• 77? // P 

2£ • V(i?Q) - — = 0, 



Rk z p*y dz 



where (Balbus 1995): 



V 



k z dz dR 



Finally, we may recombine equations ( l47j) and ( l4~8j) . separately isolating £r and ^: 
x \ o i, 2mk R - 2QRk 2 , /, • • 



m A; 2 



(4£ 



(49) 



1 ^ , /A; 2 m 2 dP\ 

+ -p {i ' Va) \t- VP -—^ ] = °' m 



R 2 k 2 dR 
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where k\ = k 2 + m 2 / 'R 2 , and 

•■ . 9> . 2mk z ■ 2QRk R k z / • , -\ 

+ (* • tu) a & - • v(izn) - — ^ + 

1 f . __ v fk R k z ^^ m 2 dP\ „ 

-^• W) (^ pp + ^&) =0 ' (51) 

Equations (l5TJj) and (l5"Tj) are the fundamental coupled equations governing the behavior of 
Lagrangian displacements. 



2.8. Plane wave limits: axisymmetry and uniform rotation 

2.8.1. Axisymmetry 

It is important to establish the axisymmetric behavior of the disturbances, since it 
represents the long time behavior of the nonaxisymmetric reponse. In particular, we have 
already noted that at large times the wavenumber ratio k R /k z is just the time-steady ratio 
of the corresponding Q gradients. By contrast, m remains fixed, so the mode becomes 
asymptotically axisymmetric as time increases and the poloidal wavenumbers grow. Thus, 
at late times all nonaxisymmetric modes behave as an axisymmetric mode whose value for 
k • v A is fixed but arbitary, while the value for k R /k z is fixed by equation (130]) . 

In the Appendix, the m — > limit of equations ( 150|) and f[5Tj) is shown to lead to the 
dispersion relation of Balbus (1995): 



or 



-Lv(R 4 n 2 ) + — {VP){Va) 
R 6 p7 



4Q 2 (k ■ v A ) 2 = 0, w 2 = co 2 - (k ■ v A ) 2 . 

(52) 



7 u; 4 + u 2 



±V(RW) + ±(VP)(Va) - 2|(fe • v A ) 2 



+ ^(k-v A ) 4 -(k-v A ) 2 



rvq 2 + — (vp)(Vo) 
pi 



(53) 



For applications to the SCZ, we are interested in the case in which the dominant balance 
of (153j) is between the first two terms, and the magnetic terms are unimportant (Goldreich 
& Schubert 1967): 

' 1 ~ ' — (VP)CDa) = (54) 

PI 
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With k R /k z = (dn/dR)(dn/dz)~\ we find that 



(55) 



an 



\~dz~) \dR~dz~~ ~dz~dR) ~ \~dz~) { >' e<p ~~fr\dl) r~d0- 



Hence, in terms of the gravitational field g = — (1/ p)(dP/dr) 



P7 7r \dz J 09 



(56) 



(57) 



The right side of this equation consists of directly observed or easily calculated quantities. 
By similar reasoning, 



^ = BMr) + , 1 = _ii(_) ._ 



(58) 



Here we have used the fact Da' = 0, since a' shares isosurfaces with Q. Combining equations 
§1D, §5}, (EH) and (PD, we obtain 



u 2 = ivn 



-2 



( —\ 2 AVL 2 + ( — \ 2 9 d(Tr 



\ dz 



d6 J 7r 2 dr 



(59) 



This dispersion relation is an interesting blend, melding a standard form for inertial/gravity 
waves in a uniformly rotating medium with k oc Vf2 for the poloidal wavenumber com- 
ponents, which obviously requires the presence of differential rotation to be sensible. The 
dynamical effects of the rotational shear are lost (DQ = 0) when k is parallel to Vf2. 



For the interesting case of Q = Q (6), the dispersion relation is 

g da r 



AVL Z sin" 9 + 



7 dr 



(60) 



In principle, these axisymmetric modes can be rotationally stabilized at equatorial or pos- 
sibly significantly higher latitudes, depending on how the adverse radial entropy gradient is 
modeled. If present, this stabilization is of some practical importance in models in which 
the Sun is convecting in surfaces of constant Q: these asymmetric modes are not unstable. 
Non-axisymmetry is more than a complication, it is crucial to the convection process itself. 



2 The R partial derivative is of course always taken with z constant, and vice-versa; the r and 9 partial 
derivatives bear a similar relationship. 
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2.8.2. Uniform rotation limit 

The dispersion relation for nonaxisymmetric plane waves in the limiting case of a uni- 
formly rotating medium may be derived from equations (150]) and (1ST]) : 



W + W 



-VP Va - 4tf) S + — ^ — VP-Vd 



4^Q 2 (k-v A ) 2 = (61) 



where ro 2 is defined in equation (I5"2l . Restricting the discussion to the nonmagnetic subscase, 
the dispersion relation becomes (Cowling 1951): 



u 2 



1 „ „ „ ,„ \ k 2 m 2 



VP Va - 4^ -§ + — — VP-Va 

P7 / Ar k z K z, yp 



0. (62) 



Notice that m introduces its own form of coupling to the entropy gradients. Moreover, it 
is possible to eliminate rotation-induced moderating influences on the growth rate by first 
setting k z = 0, and then to maximize this rate by setting k R = 0. The maximum growth 
rate is the magnitude of the Brunt- Vaisala frequency \N\, corresponding to precisely radial 
r displacements. The pure m modes are more unstable than are modes contaminated by 
poloidal wavenumber components. 

How is it possible for an element to move along a spherical radius in a rotating system 
without encountering Coriolis deflections? The answer, as noted earlier, is by striking a 
geostrophic balancqfl of the azimuthal forces: 

2pR\N\Sl£ R = -im5P (63) 

With k r = kg = 0, there are neither components of the Coriolis force in the eg or e r directions, 
nor are there unbalanced pressure gradient forces. The convective rolls occur preferentially in 
planforms of latitudinal arcs: the familiar "banana cells" often seen in laboratory experiments 
and similuations (e.g., Hart et al. 1986). Rotational forces have no effect on these rapidly 
growing linear disturbances, and therefore no effect on the convective stability criterion. 
If the squared Brunt-Vaisala frequency is negative, no amount of (uniform) rotation can 
stabilize modes with vanishing poloidal wavenumber components. 



2.8.3. Differential rotation: leading order effects 

The full problem of the evolution of three-dimensional disturbances in a two-dimensional 
background medium is a matter of some complexity, which we defer to the next section. But 



5 Pcrhaps "heliostrophic" is a more apt description for the problem at hand. 
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for applications to the SCZ, |rVlnfi| ~ 0.1, and it is appropriate to use this as a small 
parameter as a means to calculate and understand the leading order effects. If we take 
f^R = k z = as the zeroth order solution, then the polodial wavevector k p is (to all higher 
orders, in fact): 

k p = -mtVtt. (64) 

Expanding equations (150]) and (l5Tj) to linear order in the fl gradients, we obtain two very 
simple equations: 

dQ 2 • 1 dP 

in - J2— & - — — {£ • V)a = 0. (65) 
Oz 7p OR 

dQ 2 • 1 BP 
£ 2 + R—ttR - —j- (€ • V)a = 0. (66) 
oz 7p oz 

To derive equations ( 165]) and (|66|) . note that terms linear in dfl 2 /dz come from the third and 
fourth terms in equation (|50|) . and from the third term in equation (I5ip . In the former case, 
there is a cancellation of the £r terms, leaving a lone contribution from £ z . In the latter case, 
only the £r contribution is linear in the f2 gradient. In the end, only one additional term 
arises in each of equations (1651) and ( l66i) from the differential rotation, and in each case only 
the axial gradient of Q 2 enters: if the rotation is constant on clylinders (Q = Q(R)), there 
is no leading order (linear in the Q gradient) correction to the convective displacements. In 
other words, baroclinic vorticity must be present in the background to obtain a deviation from 
radial motion at linear order in the differential rotation parameter. 

It is instructive to write (|65|) and (|66|) in terms of £ r and £$, the (spherical) radial and 
colatitudinal displacements. If P ~ P(r), then to first order, fl65l) and f|66l) combine to give 

dQ 2 ■ 1 BP 

Z r + R °tZ B --Z-(t.V)a = Q, (67) 
oz 7 Or 

BVt 2 ■ 

6 - R -fa%r = 0, (68) 

Equation (1681) shows that the amplitude of will be linear in dQ 2 /dz. This means that 
in equation (1671) . all of the £g terms will be second order in the f2 2 gradient. Since we are 
working to linear order in the angular velocity gradient, our equations become 

Cr + iV 2 £r = 0, (69) 

ie - (70) 
r7 afc* 

where in equation flTOl) . we have used thermal wind balance (jlj) to substitute for Rdfl 2 /8z. 
To this order, there is no change in the behavior of the radial component of the displace- 
ment, while a poleward-increasing entropy profile produces poleward deflections of a radially 
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outward moving convective displacement (i.e., one that bears excess entropy). With |iV|i of 
order unity, the angular deflections of convective displacements from linear theory are very 
similar in magnitude to the observed departures of the iso-Q surfaces from constant 9 cones 
in the bulk of the SCZ. This, we suggest, is no coincidence: convection and its hallmark 
of constant residual entropy are both intimately associated with constant Q surfaces in the 
bulk of the SCZ (BLW), and a poleward deflection of the fluid elements is unavoidable when 
dVt/dz < 0. The interesting point, as we have earlier noted, is that it seems there must be 
an external source of vorticity in place to drive the deflections. 



3. Numerical solutions 

3.1. Representative parameters 

We next consider exact numerical solutions of equations (150]) and (l5Tj) . There are four 
important solar model parameters that need to be fixed: the two components of the Q 
gradient, and the two components of the a (entropy) gradient. The Q gradient components 
at a particular location may be read directly from the helioseismology data. The term 
da /89 then follows from the assumption of thermal wind balance. Finally, da/dr is taken 
from a published benchmark solar model (Stix 2004). Typical values for the Q gradient at 
midlatitudes near r = 0.85i? o are (Christensen-Dalsgaard & Thompson 2007): 

dlnVl dlnQ 

Thermal wind balance (i.e., vorticity conservation) then gives 

% * - 4 >< 10 ~ 6 (72) 
while a standard mixing length model (Stix 2004) gives: 

<)a -2.3 x 10- 5 (73) 



<91nr 



Thus, for n = 400nHz, 



4Q 2 ~ 2.5 x 10- n s- 2 , _9_dOj_ ^ _ g ? x 1(r i V 2 (?4) 

7r a In r 

(In the context of our model, once the gradients of Q 2 have been taken directly from the he- 
lioseismology results and the 9 derivative of a' from thermal wind balance, the r derivative of 
o' is uniquely determined from requiring counteraligned gradients of Q 2 and a'. Then, equa- 
tion ( !73|) is understood as the radial gradient of oy.) For the particular values in equations 
fT7T]) -fj74l. at latitudes less than 54°, the axisymmetric displacement solutions of equation 
( |60l are stable. 



-18 - 



3.2. Results 

We integrate equations fl50|) and ( 15 ip for a variety of different initial velocities and 
wavenumbers. Our strategy is to choose initial wavenumbers lying in the plane tangent 
to e r , with a random initial velocity direction (perpendicular to the wavenumber), and an 
initial displacement of zero. The ensuing trajectories are then followed. 

When the angular velocity is uniform and free of shear, the results are very simple: all 
trajectories rapidly become radial, regardless of their initial condition. In agreement with 
our analytic treatment, after a brief initial transient, there is no equatorial or poleward 
deflection, even with the Coriolis force. On the other hand, when an angular velocity profile 
is used that has been modeled with the helioseismology data, there are ~10% poleward 
deflections on time scales of a month or two, just as equations f )69|) and flTOj) would predict. 

Fig. 1. — Meridional slices of the Sun (r = 0.8i?©) at four different times (in units of 27r/Q) 
showing the relative orientation of the radial direction (dotted black line), entropy flux (solid 
lines), and isorotation surfaces (dotted red lines). Everywhere but at the highest latitudes, 
the entropy flux aligns itself with the isorotation surfaces within a typical mixing time. See 
text for further details. 

i = 0.1 i = l 
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Figure (1) shows this effect clearly. We work at r = O.8i? , at latitudes 15°, 30°, 45°, 
and 60°. At each latitude, the surface of the local isorotation curve is shown, edge-on in a 
meridional plane, by a dotted red line. The local radial direction is shown by the black dotted 
line. The colored solid line at each latitude is the direction of the local Eulerian perturbation 
velocity, projected in the meridional plane. Four different times are shown. In each case, the 
trajectory passes through the local isorotational surface on a time scale of a few e-foldings, 
i.e. something close to a mixing time. Near the pole — above ~ 60° — the trajectories begin 
and remain northward of the local iosortation surface, becoming more so with time. Here, the 
problem becomes close to one-dimensional, with all flow quantities predominantly a function 
of z. Our fundamental assumption that there is a functional relationship between angular 
velocity and residual entropy is then a matter of mathematical symmetry, not dynamics. 

The torques required to maintain the fluid elements in or near surfaces of constant 
angular velocity are provided by azimuthal pressure gradients, as previously noted. In no 
sense is the angular momentum of an individual fluid element conserved. Consider however 
5v^, the Eulerian perturbation of the angular velocity: 

Svj, = (R£ • V)fi = -R(Ki r + kote)/m, (75) 

where we have used equation (145]) in the last equality, and switched to spherical coordinates. 
The wavenumber component k r is very small in the bulk of the SCZ since it is proportional 
to dfl/dr and Q ~ Q(8). The term kg£g term is quadratic in the Q gradients, thus also a 
very small quantity. Therefore, there is very little angular momentum transport propagated 
by the correlated fluctuation tensor component ^ r Sv^. Indeed, noting that if equation (IMj) 
holds for the poloidal wavenumber components, then 

8V4 = -R(k r £ r + k e ie)/m = Rt£-Vtt, (76) 

which for a near vanishing Sv^ is a self-consistent indication that displacements do not stray 
from constant Q surfaces. Angular momentum could, in principle, still be directly advected 
by a nonvanishing poloidal mass flux, if one is present. 

3.3. Linear convection theory: a summary 

The linear behaviour of perturbations in a two-dimensional, stratified, differentially ro- 
tating medium is very different according to whether or not Q is dependent upon z. If O 
depends only upon R, isobaric and isochoric surfaces coincide, there is no source of baroclinic 
vorticity, and the dominant convective displacements are radial. There is neither poleward 
nor equatorial bias in the heat transport. If Q depends upon z, isobaric and isochoric sur- 
faces do not coincide, an explicit vorticity source is implicated, and outward moving, hotter 
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convective displacements deflect north or south depending upon whether f2 respectively de- 
creases or increases poleward. 

The helioseismology data in the bulk of the SCZ are in accord with this description, with 
dVl 2 jdz < and slightly upward cants to the constant f2 surfaces. The question remains, 
however, of what the cause of this gradient is. Indeed, it is a classical case of begging the 
question: the fundamental angular entropy gradient, putatively arising from the interaction 
of rotation and convection, itself requires the prexistence of an axial Q gradient — the very 
gradient the angular entropy gradient is supposed to be explaining! 

In the next section, we suggest a resolution to this problem. 



4. Vorticity generation 

4.1. Centrifugal distortion of stratified surfaces 

It has long been known that a uniformly rotating star cannot simultaneously be both 
in radiative and hydrostatic equilibrium (e.g. Schwarzschild 1958). The difficulty is that 
the rotation induces centrifugal distortions of the isothermal surfaces — polar flattening and 
equatorial bulging — which are incompatible with a vanishing radiative flux divergence. The 
dimensionless number that sets the scale for these distortions is the centrifugal parameter, 
which for a (fictional) uniformly rotating Sun takes the value: 

£ »=cS:~ L6xir (77) 

(We have used 11 = 2.5 x 10~ 6 .) This is a very small number. But we are interested in 
entropy gradients of order 

% ~ w -°- < 78 > 

and it therefore behooves us to take note of cent rifugally- induced angular gradients. 

Normally, a tiny amount of meridional circulation is enough to offset the unbalanced 
radiative heating. Consider, however, conditions at the outer edge of the radiative zone. The 
entropy equation is 



- + (..V)a 



= "(7 - 1)V-F, (79) 



where F is the radiative flux. If the divergence term on the right does not vanish at the radius 
where the entropy gradient does, there is no simple balance of the right and left sides of the 
equation. A balance could in principle be restored if a turbulent entropy flux divergence were 
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present, and this may indeed be representative of current conditions, but it is revealing to 
follow the breakdown of our simple uniform rotation model. Uneven heating would alter the 
flux until its divergence is minmized. (The time scale for this is simply the Kelvin- Helmholtz 
time; eo is not involved.) This altered flux divergence will generally be incompatible with 
coincidence of isobaric and isochoric surfaces, which uniform rotation demands. We are thus 
led to the generation of differential rotation with an axial component of the angular velocity 
gradient (baroclinicity) to maintain thermal equilibrium. Even in a more complex setting 
with a thermal entropy flux divergence present, the base of the convective zone will follow the 
the radiation flux divergence and almost certainly be more spherical than the equipotential 
surfaces. Once again baroclinic structure will be generated^. 

The helioseismology data show unambiguously that the radius at which the entropy 
gradient vanishes is symmetrically placed within a narrow band of pronounced differential 
rotation: the tachocline. Moving from this radius downward into the radiative zone, uni- 
form rotation takes over as the entropy gradient rises rapidly and meridional circulation is 
established. Moving upward into the convective zone, turbulent mixing quickly develops and 
there is an abrupt change of the character of the differential rotation: convective mixing 
leads at once to the coincidence of two important classes of surface, but not a coincidence 
that is compatible with uniform rotation. Rather, it is isorotational and residual entropy 
surfaces that now coincide. (Note the hidden but important role of nonaxisymmetry, as the 
convective modes are large m disturbances.) Turbulent convection, in this picture, does not 
generate "from scratch" the differential rotation in which it operates. Instead it reacts to, 
and reinforces, the angular velocity gradient bequeathed to it from the surface of vanishing 
entropy gradient. We suggest that this is a key ingredient to the organizational scheme of 
differential rotation in the Sun. 



4.2. Octopolar structure 

The centrifugal distortion of equipotential surfaces engendered by uniform rotation 
leaves temperature T, pressure P and density p functions of r[l + e(r)P 2 (cos 6)], where 
P2 is the usual Legendre polynomial of second order, and e is a function only of radius r (of 
order eo). The function e is determined by the solving the Poisson equation with appropriate 
boundary conditions (Schwarzschild 1958). At first glance it might be thought that the effect 
of vorticity generation and its differing iso-surfaces would be to create a distinct e for each 
structural variable: e p (r), €t(t), etc. In fact, when vorticity is generated it is impossible 



4 We acknowledge the referee M. Miesch for emphasizing this point. 
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to satisfy basic rotational equilibrium without a P4(cos#) dependence in the leading order 
nonspherical structure of these variables. 

The equation of rotational equilibrium is 

m 2 e R = -VP + V$ (80) 

p 

where Q depends upon r and 9. From the symmetry of our problem we would expect Q 2 to 
be of the form 

VI 2 = n (r) 2 + q 2 (r)P 2 (cos6) + g 4 (r)P 4 (cos 9) + ... (81) 

where the q 2 i are functions of r only. This is not necessarily an expansion in a small parame- 
ter, but for the region of interest the first two term provide a very good approximation, with 
errors in fi 2 — SIq about 10% very near the equator and less elsewher^]. The average of fi 2 
is the average of Qq and will be denoted Q 2 . Then, we may write the force balance equation 

as 

R(Q 2 -W)e R = -VP + V$' (82) 
P 

where 

$' = $ - ^-9? (83) 



M) 



Taking the divergence of equation (182 



or 



i_ a_ 

RdR 
d cot 6* 9 



P 2 (fi 2 -fi 2 ) = V- QvPJ +47rGp-2fi 2 



dr r 36 



r 2 sin 2 9(Q 2 - fi 2 ) = V- ( ~ VP ) + 4ttGp - 2fi 2 (85) 



From the form of this equation, it is evident that if fr has terms through order p(cos6*) 
in its angular expansion, then p and P will in general have angular structure through p + 2 
in the first order linear term of a small eo expansion. Therefore, uniform rotation results 
in quadrupolar deformation of the iso-surfaces of all structural variables, whereas a simple 
sin 2 9 latitudinal dependence of Q 2 results in octopolar structure. 

Consider a scenario in which the Sun is rotating uniformly. There are P 2 distortions 
of iso-surfaces in the radiative zone, and no polar deflections of warm fluid elements in the 
convective zone. Baroclinic vorticity would then be produced at the radius of vanishing 



5 Numerical fits generally are quoted for rather than SI 2 ; the latter turns out to have a simpler expansion 
near the tachocline boundary. See Gough [2007] for convenient parameterizations of 17. 
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entropy gradient. Were only P 2 structure to remain, the only self-consistent solution for 
the angular velocity would be Q = Q(r). However, dVt/dz < would now also be present, 
causing first order poleward deflections of warm convective elements. 

These results are very suggestive, and offer at least a heuristic approach to understand- 
ing the Sun's poleward decreasing angular velocity profile. If dVl 2 /dz < is maintained 
through the tachocline into the bulk of the SCZ, in which convection largely eliminates the r 
gradients of fl and a', then dVt 2 /86 must be positive. To justify these assumptions in detail, 
however, one would need to know how angular momentum is transported and deposited by 
secondary flows (Meisch et al. 2012), and what is relative importance of vorticity forcing 
by nonconservative forces versus inertial vorticity conservation. These considerations, whose 
detailed origin lies outside the the scope of the current work, regulate the the 9 dependence 
of the pressure, density, entropy, and ultimately the angular velocity. 



4.3. Solution for fi 2 

The above considerations suggest that we regard the p (for example) as a function of 
the quantity 

r p = r[l + e p (r)f p (co S 2 9)], (86) 

and similarly for P and T. The /j(cos 2 9) functions are linear combinations of P 2 and P4, in 
general distinct for i = p,P, T. (Recall that in the numerical simulations described by Miesch 
et al. (2006), P; angular structure in the entropy was included as a boundary condition at 
the base of the convection zone. By far the best solar fit included P 2 and P 4 terms in the 
angular structure of the entropy.) Expanding r p , 

p{r P ) = Po(r) + Sp(r) + ^)^-f p (co S 2 &) + ... (87) 
where p is the nonrotating solution and Sp = p(r) — Po( r )- To leading order, 

! = ^-^e^^M ,88) 

with similar results for P and T. The notation /' denotes a derivative with respect to cos 2 9. 
The /' functions are thus a linear superposition of cos 2 9 (or if more convenient sin 2 9) plus 
a constant term. 

With this development, the vorticity equation becomes 

— = 2cos9(e P f P - c ,/,)___. (89) 
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The demands of radiative equilibrium will require isothermal surfaces to be more spherical 
(less distorted) than are isochoric surfaces, whose centrifugal deformations are less critical 
to the heat flux. The pressure, being a product of p and T, would thus also be closer 
to spherically stratified than p itself. A simple mathematical approximation that allows 
progress and respects the observational fundamentals is to assume that the / stratification 
functions are the same for each variable, but that the e functions differ in magnitude from 
one variable to the next: e p > ep > We will also ignore the spatial variation of the e 
within the tachocline boundary layer. 

These considerations lead to a vorticity equation of the form 

dQ 2 A 



dz r 



cos9(sm 2 9-a) (90) 



where A (of order eoGM ) and a (of order unity) are positive constants. We have assumed 
1/r 2 gravitational field and a 1/r dependence for dlnp/dr. Signs have been chosen so that 
dfl 2 /dz is negative at non-equatorial latitudes, as indicated by the helioseismology data. 
Another way to write this equation is 

Vfl 2 A 



^r(sm 2 e-a) (91) 

where the characteristic derivative T>/T>r is taken along the path 

r 2 sin 2 9 = Constant = r 2 sin 2 9 (92) 

In this form, the equation is close to the tachocline equation (19) of BLW, which is "correct" 
(in the sense that it agrees extremely well with the data), but derived from a completely 
different point of view. The sole difference in the two formalisms is that the T>/T>r derivative 
of BLW is taken along the path 

r 2 sin 2 9 = r 2 sin 2 9 + /3r 2 (l - -) (93) 

where ft is a number of order unity. This is the characteristic path associated with the SCZ 
isorotation contours. 

How closely do the two approaches agree? If Q 2 is a given function ^o( r o ^o) & t 
radius r = r , then the solution of f l9T]) and ( l9~2"j) is 

Q 2 = ng(rg sin 2 O ) + A £ - ^ dr (94) 



or 

^2 ^2/ 2 • in \ v4sin 2 6> A Aa / 



(95) 
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Next assume that 

n 2 = + q 2 . sin 2 do ( 96 ) 

where Q 2 and Vl\ are positive constants. (Recall that this is an excellent approximation for 
the boundary of the tachocline.) Using equation (1921) and solving for the isorotation curves 
leads to 

^ = 4CV^fQ (97) 



r 2 V 1 + CF f 
where 

A = (tf-ni)/nl, C = A/(ny ), F j = (l/j)[l- (r /r)1 (98) 
Note that A is numerically the same as sin 2 #o; C is a number of order eo x 1/eo ~ unity 

As can be seen in figure (2), which shows our solution fl97|) next to the precise fit of BLW, 
the essential features of our result are basically correct. We have embedded the tachocline 
solution inside the same convective envelope solution in both cases. The principal area of 
disagreement is the bifurcation zone, on one side of which the isorotation contours break 
toward the pole and on the other side toward the equator. This is clearly the region where 
a crude approximation of the right side of equation (189]) is likely to be most inaccurate. 
Our approach to the tachocline structure makes it more clear why the BLW equation works 
so well. An uncertain approximation used by BLW was the extension of the same simple 
functional dependence between a' and Q 2 into the tachocline. But in the region of interest 
near the bifurcation point, the tachocline solution really does become a smooth extension of 
the SCZ (and the tachocline belies its name). Away from the bifurcation zone, the functional 
dependence hardly matters, as it affects only the subdominant dQ 2 /d8 term in the governing 
equation. 

The main point that one should take away from this exercise is that the combination of 
direct vorticity forcing together with a proper reckoning of the octopolar distortions of the 
stuctural variables P and p caused by the centrifugal forces of even the simplest latitude- 
dependence of the angular velocity seem to be important components of the rotational profile 
of the tachocline. 



5. Conclusions 

The current work combines two very different types of calculation, linking the linear 
dynamics of the convective zone to the origins of solar differential rotation and vorticity. We 
begin with the second part first. 



We have argued that the centrifugal distortion of equipotential surfaces combined with 
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demands of thermal equilibrium requires the cleaving of isobaric and isochoric surfaces, and 
is likely to be the underlying cause of the Sun's differential rotation. Current numerical 
simulations are not designed to capture a process in which an order eo ~ 10 ~ 5 radiative 
effect is turned into relative angular velocity gradients of order 0.1. It may well be possible, 
however, to devise other computational schemes tailored to working with this point-of-view. 

An important technical point is that the tachocline pressure and density (and therefore 
entropy) must have both P 2 (cos#) and P4(cos#) angular structure. Although this in itself is 
not a particularly new result, neither has it been widely appreciated, and we have exploited 
it in a rather novel manner. With P4 and P 2 structure in place, not only may one understand 
the simple form of the vorticity equation in the tachocline, with reasonable approximations 
one may explicitly solve the equation. An important parameter of this equation, namely the 
precise angle at which dQ 2 /dr changes sign, is probably determined by minimizing the torque 
on the radiative interior. Finally, it is interesting to note that Roxburgh (2001) calculated 
the Sun's multipole moments using models based on helioseismic inversions. He found that 
the size of the octopole term J4 was comparable to the change in the quadrupole term J 2 
when differential rotation was included. This is what one would expect if differential rotation 
were modifying P 2 and creating P 4 angular structures. 

In the first part of this paper, we have carried out a very general Lagrangian linear cal- 
culation of the fluid displacements in an arbitrary, magnetized, two-dimensional background, 
stratified in both R and z (or r and 9). Applied to conditions appropriate to the SCZ, in the 
absence of a background dVt/dz gradient, we find no poleward deflection of hot convective 
fluid elements. Such an effect is sometimes invoked to produce an angular entropy gradient 
which would in turn lead to differential rotation via thermal wind balance. Although there is 
nothing in our calculations that would prohibit the emergence of the required gradients in Vt 
and z at nonlinear order in a turbulent fluid, it is some significance that the dominant leading 
order linear response is completely different depending upon whether dfl/dz is present or 
absent. When a finite dVt/dz is a priori present, there are reinforcing convective deflections: 
a negative axial Q gradient, in particular, engenders poleward deflections of warmer fluid 
elements, which via thermal wind balance strengthen and maintain this same Q gradient. 
The presence or absence of background baroclinic vorticity is thus mirrored in the leading 
behavior of convective displacements. 

The actual generation of baroclinic vorticity may well lie outside the realm of convection 
dynamics. We propose that a vorticity source will inevitably appear at the location of a 
vanishing entropy gradient. In general, because of the effects centrifugal flattening, the 
constraints of thermal and dynamical equilibrium will force different iso-surfaces for density, 
temperature, and pressure. This must lead to an axial angular velocity gradient. By making 
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some simplifying but plausible approximations, one may calculate a time-steady solution 
of the vorticity equation for the angular velocity This explicit solution yields a tachocline 
structure that certainly resembles the observations, and bears comparison with an earlier, 
more accurate, but also more phenomeno logical, calculation (BLW). With the onset of fully 
developed convection, surfaces of constant angular velocity and residual entropy coincide, 
and the character of the isorotation contours takes on the classical conical form well-known 
from helioseismology. 

The strength of this view of the origin of solar differential is that it emerges from just a 
few rather simple and largely inevitable processes: the demands of thermal energy balance in 
a rotating system (which creates a baroclinic axial gradient of Q at the radiative/convective 
boundary), the kinematics of shear (which, via embedded wavenumbers, incorporates fl gra- 
dients into the response of the fluid displacements), and the linear dynamics of convection 
(which causes poleward displacements of entropy bearing fluid elements). While direct nu- 
merical simulation of this scenario is likely to be very challenging, it may well be possible to 
design a test-of-principle proxy system. 
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Appendix 

To recover the axisymmetric limit of equations (|50|) and (ISTl) . we begin with their equiv- 
alent forms (HU) and (135]) : 



2QR 



VP 

k R ^ R + U z ) + (£-V)<7 = 0, (99) 

PI 



Rk R 



m 



Rk z m 
+ 



m Rk 2 



i + (k ■ v A fi z 



This preprint was prepared with the A AS IATgX macros v5.2. 
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In the axisymmetric m — > limit, equation ( 199]) may be written 

A- 2 



2SIR 



rn 



k R ^ R + KCz + 



A ^ 



^ + (fe • v A y£ R 



VP 

PI 



(£-V)a = 0, 



(100) 



(101) 



In what follows, we will also require the twice differentiated form of this equation, 



OOP .. h 2 r ■■ i VP ■■ 

m k z z L J p7 



(102) 



where the notation m and iv denotes three and four time differentiations, respectively. Only 
the leading order terms in a small m expansion have been retained. 



Now, in the axisymmetric limit, equation fllOOp becomes 



R 



rn 



k R £ R + U z ) + (k • v A f- {k R £ R + k z Q - 2i ■ V(i2fi) = 



rn 



(103) 



Differentiating once 
R 



m 



{k R CS + kfi z ) + {k-v A f-(k R i R + k z i z )-R i + {k-v A fi -vn-2Z-v(Rn) = o. 

m V / L J 

(104) 

The axisymmetric dispersion relation now follows from substituting equations fllOip and 
(11021) for the 1/m terms into equation (11041) . setting £ z = —k R ^ R /k Z) and replacing all time 
derivatives by — iu. (The sign is chosen so that a positive wavenumber has a positive phase 
velocity.) After algebraic simplification, the result is 



W + W 



^V{RW) + -{VP){Vo) 



PI 



4fi 2 (fc ■ v A ) 2 = 0, w 2 = u 2 - (k ■ v A ) 2 , 

(105) 

where V is defined in equation ( I49p . This is in precise agreement with Balbus (1995). 
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Fig. 2. — Top: The C = 1, a = 0.8 tachocline solution of equation (l§Tj) along the trajectory 
characteristics of ( 1921) . Tachocline solutions are joined, at r = O.77_R , to the best fit solution 
of BLW for the SCZ and outer layers. Formal tachocline location is indicated by green line. 
Middle: The best fit global solution of BLW, in good agreement with observations. In this 
case, equation f l9T|) is solved along characteristics given by f )93|) . Bottom: GONG data, 
courtesy R. Howe. 



